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Abstract 

While navigation systems for the determination of the orbit of the Global Position System (GPS) have 
proven to be very effective, the current issues involve lowering the error in the GPS satellite ephemerides 
below their current level. In this document, the results of an orbit determination covariance assessment 
are provided. The analysis is intended to be the baseline orbit determination study comparing the benefits 
of adding laser ranging measurements from various numbers of ground stations. Results are shown for 
two starting longitude assumptions of the satellite location and for nine initial covariance cases for the 
GPS satellite state vector. 


1. Introduction 

The task of the laser ranging analysis effort is to determine the added benefits derived from solving 
for a spacecraft’s state vector when utilizing laser ranging measurements in addition to the current use of 
pseudo-range and accumulated delta range (ADR), also known as Integrated Doppler measurements. The 
methodology to perform this analysis is to perform a covariance study for the Global Positioning System 
(GPS) orbit. 

The first method will be to perform the orbit determination (OD) using pseudo-range and ADR 
measurements from the six Monitor Stations (MS) for the GPS satellite orbit through an Extended 
Kalman Filter (EKF) analysis. The second method will be to add laser ranging measurements from the six 
MS sites along with various amounts of additional sites. Measurements from both methods will be used to 
form estimates for the satellite’s state vector, which will be propagated until new measurements are 
available. Finally, comparisons will be made between the performance of the current system and the 
modified system. 


2. EKF Description 

The purpose of an Extended Kalman Filter (EKF) is to estimate the states of a non-linear system. The 
EKF is an extension of the standard Kalman Filter, where its purpose is to estimate the states of a linear 
system. The derivation of the EKF is based on linearizing the non-linear system using the Kalman Filter 
estimate as the nominal state trajectory. The non-linear system is linearized around the Kalman Filter 
estimate and the Kalman Filter estimate is based on the linearized system (ref. 1). 

The method of the EKF used for these simulations is the discrete time system/discrete time 
measurement EKF. This is the most appropriate method to simulate the EKF, as performing continuous 
time dynamics in a computer requires an extremely large amount of memory and processor power to be 
performed efficiently. Also, it is important to note up front that typically, there are multiple runs 
performed for each scenario. This is due to the fact that the equations of the EKF dictate that the real 
noise parameters, instead of the covariance of the noise (as in linear Kalman filter problems), be used to 
form new estimates of the state estimates. The following equations dictate the nature of the EKF process 
(ref. 1). 
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1) System and Measurement equations 


x k = fk-x{ x k-X,U k -X,Wk-\ 
y k = h k (x k ,v k ) 
w k ->N(0,Q k ) 
v k — > N{0, R k ) 

Where: 

• fk is the state equation at time k 

• x k is the state vector at time k 

• u k is the input vector at time k 

• w k is the state noise vector at time k 

• h k is the measurement equation at time k 

• y k is the measurement vector at time k 

• v k is the measurement noise vector at time k 

• Q k is the state noise covariance at time k 

• R k is the measurement noise covariance at time k 

2) Filter Initial Conditions 

*o = £ ( x o) 


Po +=E 


(*o - x 0 to - X 


Where: 

• e( ) is the expectation operator 

• Xq is the a posteriori state estimate at time 0 

• Pq is the a posteriori covariance estimate at time 0 


3) Filter Loop 

a) Compute the following partial derivative matrices 


F k - 1 


dfk-X , 
dx x k- x 


L k -x 


dfk-x , 
dw x k - 1 
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Where: 


df 

— | is the partial derivative operator of function /by vector x, evaluated by vector x, 

fix 1 

F k is the partial derivative matrix of the state equation by the state vector at time k 
L k is the partial derivative matrix of the state equation by the state noise vector at time k 


b) Perform the time update of the state estimate and estimation error covariance 

P k = F k-\ P k-\ F k-\ + L k-\Qk-\L T k-\ 

*k = fk-l{x£-l’ U k-l’°) 


( 9 ) 

( 10 ) 


Where: 

• P k is the a priori covariance estimate at time k 

• x k is the a priori state estimate at time k 


c) Compute the following partial derivative matrices 

fix X k 


HI) 


M k = 


dhk_ 

dv 



( 12 ) 


Where: 

• H k is the partial derivative matrix of the measurement equation by the state vector at time k 

• M k is the partial derivative matrix of the measurement equation by the measurement noise vector 
at time k 

d) Perform the measurement update of the state estimate and estimation error covariance 

K k = P k ~HT (H k P k ~Hl + M k R k M\ )~ l (13) 

x k= x k + K k [yk - h [x k ; 0 )] ( 14 ) 

P k ={l ~ K k H k )P k (15) 


Where: 

• K k is the Kalman gain matrix at time k 

• x£ is the a posteriori state estimate at time k 

• Py is the a posteriori covariance estimate at time k 
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3. State Dynamics Description 

The orbit of a GPS satellite can be viewed as a simplified two-body problem. That is the assumption 
made for this analysis. Therefore, the differential equation that solves the two-body problem is as follows, 
in Eq. (16) (ref. 2). 


r + ^ r = ad 


(16) 


Where: 

• r is the position matrix 

• r is the 2 nd time derivative of the position matrix 

• p is the Earth’s Gravitational constant 

• R is the magnitude of the position matrix 

• is the orbital perturbation 

It is important to note initially that the OD analysis solves for more than just the position matrix. The 
purpose of the OD analysis is to solve for position, velocity, clock bias, and frequency bias. Therefore, the 
state equation that governs the OD analysis is an extension of the two-body problem. Keep in 
consideration that the orbital perturbation is assumed to be zero for this analysis. The state is defined in 
Eq. (17), while the state equation is given as follows in Eq. (18). 


x = 


r 

v 

Ctbias 

Ctbias 


(17) 


V 


' 0 ~ 

-p 


0 

R 3 

x + 

0 


1 

0 


.1 // 


(18) 


Where: 

• r is the position matrix 

• v is the velocity matrix 

• c is the speed of light in a vacuum 

• tbias is the clock difference between the satellite and the ground stations 

• /is the GPS L 1 frequency 

• w is the state noise 

Equation (19) shows the fully expanded form of the discrete state equation. Equation (20) converts 
the state equation from Eq. (18) into a discrete time state equation, needed for the discrete time EKF, in 
Earth-Centered Fixed coordinates. 
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x k, 1 


r k,x 

x k,2 


r k,y 

x k,T> 


r k,z 

x kA 


X k,x 

x k, 5 


Vk,y 

x k, 6 


Vk,z 

x k,l 


Ctbias 

x k, 8 _ 


Ctbias _ 


( 19 ) 


x k+\ 


cos(c))) 

sin(c))) 

0 

Acos(i|)) 

Asin(cj)) 

0 

- sin(c))) 

cos(i|)) 

0 

— A sin(c()) 

Acos(i))) 

0 

0 

0 

1 

0 

0 

A 

-^T cos W 

-^T sin W 

0 

COs(())) 

sin(i|)) 

0 

^T sin W 

-^y C0S W 

0 

- sin(c))) 

cos(i|)) 

0 

0 

0 

_ JL 

0 

0 

1 

R 3 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

1 

0 


\Xk + 


0 

0 

0 

0 

0 

0 

1 

1 // 


Wk 


( 20 ) 


Where: 

• ()) is the Earth rotation rate in radians/second 

• A is the discrete time step 


4. Measurement Description 

There are three different measurement types that are utilized within the trade space of this analysis. 
The first measurement is the pseudo-range (PR) measurement, which is utilized at the MS locations. The 
equation for the pseudo-range measurement is given in Eq. (21) (ref. 3). 

PR = -*2 ) 2 + (jl -J 2) 2 +(*i — z 2 ) 2 + ct bias + v PR (21) 


Where: 

• PR is the pseudo-range measurement 

• (x] ,y\,z\) is the position of the transmitter (or receiver) 

• ( x 2 > }’ 2 , z 2 ) is the position of the receiver (or transmitter) 

• v pr is the noise term in the pseudo-range measurement 
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The second measurement type that is utilized at the six MS locations is the Accumulated Delta Range 
(ADR) measurement. This measurement also is called carrier phase. It is the integral of the range-rate 
measurement used with instantaneous Doppler shift. Equation (22) provides the mathematical description 
of the ADR measurement (ref. 3). 


adr=£ 


k = o 


W -*2 ,a-X v *u- -v* 2 t)+(tu -y 2 ,k\m,k -vyy^+Wk ~ z 2 ,k\ vz u ~ vz i,k\ 

■J(xi,k ~*2 ,k ) 2 +{y\,k —yi,k ) 2 +( z u ~ z 2,k ) 2 


+ fat bias T V ADR 


( 22 ) 


Where: 

• ADR is the accumulated delta range measurement 

• (yx \ , vyi , vz| ) is the velocity of the transmitter (or receiver) 

• (vx 2 , vy 2 , VZ 2 ) is the velocity of the receiver (or transmitter) 

• v^dr is the noise term in the accumulated delta range measurement 

The final measurement type that is utilized only in the modified systems, but at all ground station 
sites, is the laser ranging (LR) measurement. This can be thought of as the equivalent of a two-way 
radiometric signal, in terms of the equation governing the measurement. Equation (23) provides the 
mathematical description of the LR measurement (ref. 3). 


LR = 


~x 2 ) 2 + (yj - y 2 Y + (d - z 2 Y + 


vlr 


(23) 


Where: 

• LR is the laser ranging measurement 

• v LR is the noise term in the laser ranging measurement 

Table 1 provides the standard deviation of the noise terms that are assumed for the three measurement 
equations provided in Eqs. (21) through (23). 


TABLE 1.— STANDARD DEVIATION OF 
MEASUREMENT NOISE TERMS 


Noise term 

a 

V PR 

2 m 

v ADR 

5 mm/s 

V LR 

lm 
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5. Station Locations 


The first of the two proposed systems that are analyzed is the current system, which utilizes pseudo- 
range and ADR measurements from the six MS locations. The second proposed system is the modified 
system, which utilizes all measurements for the current system, plus laser ranging measurements from the 
six MS locations along with measurements from various amounts of additional ground stations. The 
locations of the additional ground stations are from a listing of laser ranging sites from the International 
Laser Ranging Service (ref. 4). The six MS have their locations listed in table 2, below. Figure 1 
illustrates the locations of the six MS sites, shown in red dots. 


TABLE 2.— MS LOCATIONS 


Station 

Longitude 

Latitude 


(°E) 

(°N) 

Colorado Springs 

-104.7167° 

38.8167° 

Hawaii 

-158.2519° 

21.5653° 

Cape Canaveral 

-80.5333° 

28.4667° 

Ascension Island 

-14.3700° 

-7.9500° 

Diego Garcia 

72.4000° 

-6.5667° 

Kwajalein 

167.7319° 

8.7189° 



Figure 1 . — Current system station map. 
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There are five modified system concepts. Modified System 1 (MSysl) utilizes all measurements from 
the current system plus laser ranging measurements from the six MS locations and two additional 
locations. Table 3 lists the locations for the two additional ground sites used for laser ranging 
measurements only for MSysl. Figure 2 illustrates the locations of the six MS and the two additional 
ground stations. The six MS are shown in red dots, while the two additional ground stations are shown in 
blue dots for the MSys 1 . 


TABLE 3.— ADDITIONAL LASER RANGING 


STATION LOCATIONS— MSysl 


Station 

Longitude 

(°E) 

Latitude 

(°N) 

GSFC 

-76.8760° 

39.0044° 

McDonald Observatory 

-104.0217° 

30.6717° 



Figure 2. — MSysl station map. 
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Modified System 2 (MSys2) utilizes all measurements from MSysl plus laser ranging measurements 
from two additional locations. Table 4 lists the locations for the two additional ground sites used for laser 
ranging measurements only for MSys2. Figure 3 illustrates all of the locations of the ground sites for 
MSys2. The eight red dots represent stations from MSysl, while the two additional ground stations are 
shown in blue dots. 


TABLE 4.— ADDITIONAL LASER RANGING 
STATION LOCATIONS— MSys2 


Station 

Longitude 

(°E) 

Latitude 

(°N) 

Hartebeesthoek 

27.6861° 

-27.6861° 

Wuhan 

114.4897° 

30.5157° 



Figure 3. — MSys2 station map. 
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Modified System 3 (MSys3) utilizes all measurements from MSys2 plus laser ranging measurements 
from four additional locations. Table 5 lists the locations for the four additional ground sites used for laser 
ranging measurements only for MSys3. Figure 4 illustrates all of the locations of the ground sites for 
MSys3. The 10 red dots represent stations from MSys2, while the four additional ground stations are 
shown in blue dots. 


TABLE 5.— ADDITIONAL LASER RANGING 
STATION LOCATIONS— MSys3 


Station 

Longitude 

(°E) 

Latitude 

(°N) 

Concepcion 

-73.0253° 

-36.8430° 

Mt. Stromlo 

149.0099° 

-35.3161° 

Komsomolsk 

136.7438° 

50.6946° 

Metsahovi 

24.3946° 

60.2172° 



Figure 4. — MSys3 station map. 
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Modified System 4 (MSys4) utilizes all measurements from MSys3 plus laser ranging measurements 
from four additional locations. Table 6 lists the locations for the four additional ground sites used for laser 
ranging measurements only for MSys4. Figure 5 illustrates all of the locations of the ground sites for 
MSys4. The 14 red dots represent stations from MSys3, while the four additional ground stations are 
shown in blue dots. 


TABLE 6.— ADDITIONAL LASER RANGING 


STATION LOCATIONS— MSys4 


Station 

Longitude 

(°E) 

Latitude 

(°N) 

Tahiti 

-149.6063° 

-17.5768° 

Yarragadee 

115.3467° 

-29.0464° 

San Fernando 

-6.2055° 

36.4650° 

Urumqi 

87.7100° 

43.8100° 



Figure 5. — MSys4 station map. 
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Modified System 5 (MSys5) utilizes all measurements from MSys4 plus laser ranging measurements 
from four additional locations. Table 7 lists the locations for the four additional ground sites used for laser 
ranging measurements only for MSys5. Figure 6 illustrates all of the locations of the ground sites for 
MSys5. The 18 red dots represent stations from MSys4, while the four additional ground stations are 
shown in blue dots. 


TABLE 7.— ADDITIONAL LASER RANGING 
STATION LOCATIONS— MSys5 


Station 

Longitude 

(°E) 

Latitude 

(°N) 

Wrightwood 

-117.6830° 

34.3820° 

Riyadh 

46.4004° 

24.9102° 

Arequipa 

-71.4930° 

-16.4657° 

Koganei 

139.4879° 

35.7102° 



Figure 6. — Modified system 5 station map. 
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6. Methodology Description 

The analysis is performed using the discrete time/discrete measurement EKF procedure described 
previously. The state is propagated at a rate of 1 Hz. The total simulation is set to run for 1 day, or 86400 
seconds. It should be noted however, that measurements are not taken on this same one second time 
period. Pseudo-range and laser ranging measurements are taken the instant that the satellite can see the 
ground station, and again every 60 seconds after the initial measurement until the ground station is no 
longer visible. 

For ADR measurements, the integration of the delta-range measurements begins the first second that 
the ground station is visible to the satellite. The measurement is completed once the ground station has 
been in view for 60 seconds. When the satellite loses visibility to the ground station, the integration 
procedure is reset and the satellite must start over on the next visible pass. 

Note that in order for the ground station to be visible to the satellite, the ground station must see the 
satellite with at least a minimum elevation angle of 10°. This rule applies to all three measurement types 
that have been described. 

The satellite is modeled in the GPS 12 hour orbit at an inclination of 55°. The orbit is assumed to start 
on the plane of the Earth’s Equator. However, due to the nature that the ground stations are oriented on 
the surface, the satellite is modeled on starting longitudes of 0° and 90°E. Figures 7 and 8 show the points 
over the Earth surface in which the orbit for these two starting longitudes pass. Note that the point shown 
in red is the first point of the orbit in the simulation. 


Satellite Ground Path - Starting Longitude 0° 



Figure 7. — Satellite ground path — starting longitude 0°. 
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Satellite Ground Path - Starting Longitude 90 



Figure 8. — Satellite ground path — starting longitude 90°. 


Initial errors on the order of 1 km are added to each Cartesian dimension along with 1 m per second 
velocity errors in each Cartesian dimension. Clock bias and frequency bias states also begin with a 1 km 
and 1 m per second error, respectively. The other condition that is varied for the simulation is the initial 
covariance estimate. Initial conditions for the covariance estimate are formed into nine cases, as listed in 
table 8. 


TABLE 8.— INITIAL COVARIANCE CONDITIONS 


Case 

Position 

Velocity 

1 

0.1 m 2 

0.1 (m/s) 2 

2 

0.1 m 2 

1 (m/s) 2 

3 

0.1 m 2 

10 (m/s) 2 

4 

1 m 2 

0.1 (m/s) 2 

5 

lm 2 

1 (m/s) 2 

6 

1 m 2 

10 (m/s) 2 

7 

10 m 2 

0.1 (m/s) 2 

8 

10 m 2 

1 (m/s) 2 

9 

10 m 2 

10 (m/s) 2 


Finally, for the purpose of this analysis, there will be 10 runs performed for every starting 
longitude/initial covariance case simulation. Performance along the multiple runs will be combined to 
attain overall performance. 
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7. Results 


Results for the analysis of the comparison of the six systems are shown through the use of the 
covariance estimate of the systems. Even though the variance of the noise terms is constant throughout 
the simulations, multiple noise profiles are needed to analyze the performance of an EKF simulation. This 
is due to the fact that the equations of the EKF dictate that the real noise parameters, instead of the 
covariance of the noise, are used to form new estimates of the state. This is one difference between linear 
Kalman filtering and Extended Kalman Filtering. The metric to compare the performance of the six 
systems is range covariance, which is based on the final covariance estimate at the end of the simulation. 

Equation 24 is used to compare the systems for the range covariance statistic for the covariance 
estimate. 


COV, 


RANGE ~ 


nlong * nrun r~f rrl 
* 1=1 u = l 


nlong nrun 

X Xr + P 2 (2,2),kf ,i,H + P 2 (33),kfJM 


(24) 


Where: 

• COV range is the covariance range error for the covariance estimate 

• ’ p (2,2 w,i,ii , P(3,3),kf,i,ii ) are thc covariance terms for the individual Cartesian 

dimension position terms at final time kf for run ii at longitude case i 

• nlong is the number of starting longitude scenarios 

• nrun is the number of noise profile runs 

Equation (24) is used to derive performance characteristics for the six systems that are examined. 
Results are shown in a vertical bar graph with the system on the x-axis and the range covariance statistic 
on the v-axis. These bar graphs are produced a total of nine times, to account for each of the nine initial 
covariance cases. 

Figure 9 shows the performance of the six systems under the first case of the initial covariance 
estimates for range covariance performance. This is when the initial covariance parameters are on the 
order of 0.1 m 2 and 0.1 (m/s) 2 . 


Range Covariance - Initial Covariance Case 1 



CS MSysl MSys2 MSys3 MSys4 MSys5 
System 

Figure 9. — Range covariance — initial covariance Case 1. 
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Figure 10 shows the performance of the six systems under the second case of the initial covariance 
estimates for range covariance performance. This is when the initial covariance parameters are on the 
order of 0.1 m 2 and 1 (m/s) 2 . 

Figure 1 1 shows the performance of the six systems under the third case of the initial covariance 
estimates for range covariance performance. This is when the initial covariance parameters are on the 
order of 0.1 m 2 and 10 (m/s) 2 . 


Range Covariance - Initial Covariance Case 2 
0.35 | , , , , , 


0.3 - 



CS MSysl MSys2 MSys3 MSys4 MSys5 
System 

Figure 10. — Range covariance — initial covariance Case 2. 


Range Covariance - Initial Covariance Case 3 
1 .4 1 1 1 1 1 

1.2 - 



CS MSysl MSys2 MSys3 MSys4 MSys5 
System 

Figure 11. — Range covariance — initial covariance Case 3. 
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Figure 12 shows the performance of the six systems under the fourth case of the initial covariance 
estimates for range covariance performance. This is when the initial covariance parameters are on the 
order of 1 m 2 and 0.1 (m/s) 2 . 

Figure 13 shows the performance of the six systems under the fifth case of the initial covariance 
estimates for range covariance performance. This is when the initial covariance parameters are on the 
order of 1 m 2 and 1 (m/s) 2 . 


Ranqe Covariance - Initial Covariance Case 4 
0.7, , , , , , 



CS MSysl MSys2 MSys3 MSys4 MSys5 
System 

Figure 12. — Range covariance — initial covariance Case 4. 


Range Covariance - Initial Covariance Case 5 



CS MSysl MSys2 MSys3 MSys4 MSys5 
System 

Figure 13. — Range covariance — initial covariance Case 5. 
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Figure 14 shows the performance of the six systems under the sixth case of the initial covariance 
estimates for range covariance performance. This is when the initial covariance parameters are on the 
order of 1 m 2 and 10 (m/s) 2 . 

Figure 15 shows the performance of the six systems under the seventh case of the initial covariance 
estimates for range covariance performance. This is when the initial covariance parameters are on the 
order of 10 m 2 and 0. 1 (m/s) 2 . 


Range Covariance - Initial Covariance Case 6 



System 

Figure 14. — Range covariance — initial covariance Case 6. 


Range Covariance - Initial Covariance Case 7 
0.8 1 1 1 1 1 



CS MSysl MSys2 MSys3 MSys4 MSys5 
System 


Figure 15. — Range covariance — initial covariance Case 7. 
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Figure 16 shows the performance of the six systems under the eighth case of the initial covariance 
estimates for range covariance performance. This is when the initial covariance parameters are on the 
order of 10 m 2 and 1 (m/s) 2 . 

Figure 17 shows the performance of the six systems under the ninth case of the initial covariance 
estimates for range covariance performance. This is when the initial covariance parameters are on the 
order of 10 m 2 and 10 (m/s) 2 . 


Range Covariance - Initial Covariance Case 8 



System 

Figure 16. — Range covariance — initial covariance Case 8. 

Range Covariance - Initial Covariance Case 9 
0.8 1 1 1 1 1 1 



CS MSysl MSys2 MSys3 MSys4 MSys5 
System 


Figure 17. — Range covariance — initial covariance Case 9. 
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Table 9 lists the final range covariance statistics broken out by longitude versus initial covariance 
case. The range covariance statistics are the average of the final range covariance statistics for each run. 


TABLE 9.— RANGE COVARIANCE STATISTICS BY LONGITUDE 



System 

Case 1 

Case 2 

Case 3 

Case 4 

Case 5 

Case 6 

Case 7 

Case 8 

Case 9 

0° 

CS 

0.291608 

0.45046 

1.510936 

0.66386 

2.436703 

0.84973 

0.883964 

0.884283 

1.130536 


0.056013 

0.104778 

0.051748 

0.065541 

0.015729 

0.186435 

0.081631 

0.198967 

0.009543 


0.03981 

0.0774 

0.122318 

0.044429 

0.24755 

0.115746 

0.045028 

0.045007 

0.172855 


0.030392 

0.030389 

0.030389 

0.03339 

0.033389 

0.033389 

0.033764 

0.033764 

0.033764 

[3EUSH 

002429 

0024288 

0.024288 

0025982 

0.025981 

0 025979 

0.026182 

0.026181 

0.02618 


0.019915 

0.019914 

0.019878 

0.020988 

0.020987 

0.020985 

0.02111 

0.02111 

0.021109 

90 °E 

Ics 

0.398027 

0.076659 

0.534994 

0.545768 

0.530322 

0.69956 

0.628294 

0.61492 

0.298915 


0.049903 

0 129587 

0.009439 

0048682 

0.090065 

0030502 

0.169729 

0.044325 

0.054874 


0.046169 

0.041485 

0.012127 

0.055961 

0.012876 

0.05139 

0.057204 

0.057459 

0.056998 


0.028958 

0028939 

0.026205 

0.015941 

0.030104 

0.031659 

0.031912 

0.015938 

0.032032 

cbeb 

0.022156 

0.021378 

0.019173 

0.01187 

0.021379 

0.019558 

0.023869 

0.01186 

0.02387 


0.017915 

0.016662 

0.011083 

0.011044 

0.013626 

0.016653 

0.018862 

0.010178 

0.012779 


Profiles for the range covariance parameters are shown in Appendix A. It is important to notice how 
the starting longitude and initial starting covariance affect the performance for each of the systems. It is 
easy to visualize from the graphics that the initial covariance cases in which the velocity terms have a 
lower/equal initial covariance than the position terms perform more consistently, compared to when the 
velocity terms have larger initial covariance terms than the position terms. 

It should also be noted at this time about an issue that arose during the simulations. Other starting 
longitudes of 180°E and 90°W were run, but there was difficulty with the numerical stability of the 
simulations. The difficulty was repeatability of the simulation results on different computers. 

Duplications of the simulations were run on various computers with common noise profiles and other 
parameters, with differing results. Issues that were seen included non-symmetric covariance matrices and 
differences in the inversions of ill-conditioned matrices. Possible resolution of this issue will be discussed 
in the Future Work section. 


8. Conclusions 

Results for the comparison between the current tracking system utilizing pseudo-range and ADR 
measurements from the six MS locations with the five modified systems including laser ranging 
measurements from the same ground sites plus additional sites have been provided. All of the results 
show benefits of having laser ranging measurements used to solve for the satellite’s position component 
of the state vector. The results show an initial dependency on the initial longitude of the orbit. A second 
parameter that has been shown to affect performance is the initial covariance for the system. However, for 
both of these parameters, the final covariance is not strongly affected. 

The parameter that does strongly affect the final time covariance is the number of laser ranging 
ground stations used. As seen in the results in figures 9 through 17, it is typical that with an increase in 
the number of laser ranging ground stations, the final time range covariance statistic decreases. It is 
important to note that as more and more ground stations are added to the scenario, that the final time 
covariance does not keep decreasing more and more. This means that after a certain point, there exists a 
value of the number of laser ranging ground stations for which there is little decrease in the final time 
range covariance. The plots that are provided in Appendix A illustrate this effect, even when the y-axis is 
plotted in log scale. 

The initial covariance of the state is an estimate for how well the state is understood. Typically, when 
the state’s covariance is larger, then more emphasis is placed on the measurements when producing the 
EKF Kalman gain. However, when dealing with an OD type of analysis similar to the one performed, 
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where and when there are few measurements available to the receiver, then covariance parameters can 
increase tremendously. Also, the covariance of the measurement noise is an important parameter for how 
the covariance is propagated. 

Therefore, in order to better estimate the performance of the system, there needs to be more tracking 
measurements available to be used in the EKF, especially in the beginning time periods of the simulation. 
This can be provided in two forms. First, there can be more tracking sites used to provide measurements. 
Second, measurements can be taken more frequently. Modifying these options was not analyzed in this 
initial study. The first option would mean removing ground stations at some point in time in the 
simulation, where this simulation assumed that the ground station would be present throughout the 
entirety of the simulation. Also, variations in measurement noise covariance were not considered for this 
study. 

Results shown from this study include the fact that there are differences in performance between the 
current system and the modified systems including laser ranging measurements. Performance is 
dependent on the location of the ground stations and how those stations are viewed by the satellite. 
Therefore, if the additional ground stations for the modified systems were picked differently, then the 
results would vary. Therefore, elevation angle that is seen from the ground to the satellite affects 
performance. It is believed that if constraints (such as range and/or speed) are placed within the EKF, 
performance of the two systems may be better modeled. 

The question that needs to be answered from the analysis is if laser ranging measurements are needed, 
and if so, then how many stations are necessary. First off, the answer depends on the requirements. In 
stating this, it is acknowledged that the requirements will state a certain covariance level needed. This 
should be the ultimate answer as to what system to use. However, this analysis has shown that laser 
ranging measurements are beneficial and reduce the steady state system performance. Typically, the more 
stations that are added to the scenario, the lower the steady state system performance will be. However, 
there is little change in system performance when going from MSys3 to MSys4 and then from MSys4 to 
MSys5. Therefore, given the orientation of the ground stations as such as in this report, it appears that 
MSys3 would give the most benefit. Therefore, it is believed that if laser ranging measurements would in 
the future be taken into account for doing OD analysis on a GPS orbit, then the recommendation would be 
have measurements taken from the six MS stations with measurements from an additional eight ground 
stations around the world. 


9. Future Work 

Future work in the area of OD determination for a GPS orbit comparing the current system of pseudo- 
range and ADR measurements with modified systems using varying numbers of laser ranging ground 
stations, pseudo-range and ADR measurements can be expanded in multiple methods. The following list 
is not a complete list on what can be performed to enhance the analysis provided here, but is what the 
author believes to be an important next step. 

• Model scenario with a form of square root/U-D filtering to improve numerical stability issues that 
arise from using finite numerical processing 

• Model scenario with a 2 nd order EKF 

• Perform parametric analysis with variations in how often measurements are taken 

• Perform parametric analysis with variations in measurement covariance 

• Perform further parametric analysis on starting latitude and longitude 

• Perform smoothing on all past state estimates for the purpose of obtaining a better estimate of 
new state estimates 

• Perform smoothing on all state estimates for the purpose of obtaining a better estimate all state 
estimates 

• Include orbital perturbations to make scenario more realistic 
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Appendix A — Range Covariance Results 

A.l. Initial Covariance Case 1 Results 

Figures A. 1 . 1 and A. 1 .2 show the range covariance results for the 0° and 90°E starting longitude 
conditions, respectively, when the initial covariance is operating in Case 1 parameters (Cartesian position 
components = 0.1 m 2 , Cartesian velocity components = 0.1 (m/s) 2 ). Blue lines represent range covariance 
results for the current system based on pseudo-range and ADR measurements from the 6 MS ground 
stations. Red lines represent range covariance results for MSysl, combining the current system with laser 
ranging measurements from the 6 MS stations and 2 additional ground stations. Similarly, black lines 
represent range covariance results for MSys2; green lines represent range covariance results for MSys3; 
cyan lines represent range covariance results for MSys4; yellow lines represent range covariance results 
for MSys5. Note that there are 10 individual lines for each color, derived from the 10 noise profile runs. 



Figure A. 1.1. — Starting longitude 0° results. 



Figure A. 1.2. — Starting longitude 90° results. 
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A.2. Initial Covariance Case 2 Results 


Figures A.2.1 and A.2.2 show the range covariance results for the 0° and 90°E starting longitude 
conditions, respectively, when the initial covariance is operating in Case 2 parameters (Cartesian position 
components = 0.1 m 2 , Cartesian velocity components = 1 (m/s) 2 ). Blue lines represent range covariance 
results for the current system based on pseudo-range and ADR measurements from the 6 MS ground 
stations. Red lines represent range covariance results for MSysl, combining the current system with laser 
ranging measurements from the 6 MS stations and 2 additional ground stations. Similarly, black lines 
represent range covariance results for MSys2; green lines represent range covariance results for MSys3; 
cyan lines represent range covariance results for MSys4; yellow lines represent range covariance results 
for MSys5. Note that there are 10 individual lines for each color, derived from the 10 noise profile runs. 



Figure A.2.1. — Starting longitude 0° results. 



Figure A.2.2. — Starting longitude 90° results. 
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A.3. Initial Covariance Case 3 Results 


Figures A.3.1 and A.3.2 show the range covariance results for the 0° and 90°E starting longitude 
conditions, respectively, when the initial covariance is operating in Case 3 parameters (Cartesian position 
components = 0.1 m 2 , Cartesian velocity components =10 (m/s) 2 ). Blue lines represent range covariance 
results for the current system based on pseudo-range and ADR measurements from the 6 MS ground 
stations. Red lines represent range covariance results for MSysl, combining the current system with laser 
ranging measurements from the 6 MS stations and 2 additional ground stations. Similarly, black lines 
represent range covariance results for MSys2; green lines represent range covariance results for MSys3; 
cyan lines represent range covariance results for MSys4; yellow lines represent range covariance results 
for MSys5. Note that there are 10 individual lines for each color, derived from the 10 noise profile runs. 



Figure A.3.1. — Starting longitude 0° results. 



Figure A.3.2. — Starting longitude 90° results. 
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A.4. Initial Covariance Case 4 Results 


Figures A.4.1 and A.4.2 show the range covariance results for the 0° and 90°E starting longitude 
conditions, respectively, when the initial covariance is operating in Case 4 parameters (Cartesian position 
components = 1 m 2 , Cartesian velocity components = 0.1 (m/s) 2 ). Blue lines represent range covariance 
results for the current system based on pseudo-range and ADR measurements from the 6 MS ground 
stations. Red lines represent range covariance results for MSysl, combining the current system with laser 
ranging measurements from the 6 MS stations and 2 additional ground stations. Similarly, black lines 
represent range covariance results for MSys2; green lines represent range covariance results for MSys3; 
cyan lines represent range covariance results for MSys4; yellow lines represent range covariance results 
for MSys5. Note that there are 10 individual lines for each color, derived from the 10 noise profile runs. 



Figure A.4.1. — Starting longitude 0° results. 



Figure A.4.2. — Starting longitude 90° results. 
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A.5. Initial Covariance Case 5 Results 


Figures A.5.1 and A.5.2 show the range covariance results for the 0° and 90°E starting longitude 
conditions, respectively, when the initial covariance is operating in Case 5 parameters (Cartesian position 
components = 1 m 2 , Cartesian velocity components = 1 (m/s) 2 ). Blue lines represent range covariance 
results for the current system based on pseudo-range and ADR measurements from the 6 MS ground 
stations. Red lines represent range covariance results for MSysl, combining the current system with laser 
ranging measurements from the 6 MS stations and 2 additional ground stations. Similarly, black lines 
represent range covariance results for MSys2; green lines represent range covariance results for MSys3; 
cyan lines represent range covariance results for MSys4; yellow lines represent range covariance results 
for MSys5. Note that there are 10 individual lines for each color, derived from the 10 noise profile runs. 



Figure A.5.1. — Starting longitude 0° results. 



Figure A.5.2. — Starting longitude 90° results. 
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A.6. Initial Covariance Case 6 Results 


Figures A.6.1 and A.6.2 show the range covariance results for the 0° and 90°E starting longitude 
conditions, respectively, when the initial covariance is operating in Case 6 parameters (Cartesian position 
components = 1 m 2 , Cartesian velocity components = 10 (m/s) 2 ). Blue lines represent range covariance 
results for the current system based on pseudo-range and ADR measurements from the 6 MS ground 
stations. Red lines represent range covariance results for MSysl, combining the current system with laser 
ranging measurements from the 6 MS stations and 2 additional ground stations. Similarly, black lines 
represent range covariance results for MSys2; green lines represent range covariance results for MSys3; 
cyan lines represent range covariance results for MSys4; yellow lines represent range covariance results 
for MSys5. Note that there are 10 individual lines for each color, derived from the 10 noise profile runs. 



Figure A.6.1. — Starting longitude 0° results. 



Figure A.6.2. — Starting longitude 90° results. 
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A.7. Initial Covariance Case 7 Results 


Figures A.7.1 and A.7.2 show the range covariance results for the 0° and 90°E starting longitude 
conditions, respectively, when the initial covariance is operating in Case 7 parameters (Cartesian position 
components =10 nr, Cartesian velocity components = 0.1 (m/s) 2 ). Blue lines represent range covariance 
results for the current system based on pseudo-range and ADR measurements from the 6 MS ground 
stations. Red lines represent range covariance results for MSysl, combining the current system with laser 
ranging measurements from the 6 MS stations and 2 additional ground stations. Similarly, black lines 
represent range covariance results for MSys2; green lines represent range covariance results for MSys3; 
cyan lines represent range covariance results for MSys4; yellow lines represent range covariance results 
for MSys5. Note that there are 10 individual lines for each color, derived from the 10 noise profile runs. 



Figure A.7.1. — Starting longitude 0° results. 



Figure A.7.2. — Starting longitude 90° results. 
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A.8. Initial Covariance Case 8 Results 


Figures A.8.1 and A.8.2 show the range covariance results for the 0° and 90°E starting longitude 
conditions, respectively, when the initial covariance is operating in Case 8 parameters (Cartesian position 
components =10 nr, Cartesian velocity components = 1 (m/s) 2 ). Blue lines represent range covariance 
results for the current system based on pseudo-range and ADR measurements from the 6 MS ground 
stations. Red lines represent range covariance results for MSysl, combining the current system with laser 
ranging measurements from the 6 MS stations and 2 additional ground stations. Similarly, black lines 
represent range covariance results for MSys2; green lines represent range covariance results for MSys3; 
cyan lines represent range covariance results for MSys4; yellow lines represent range covariance results 
for MSys5. Note that there are 10 individual lines for each color, derived from the 10 noise profile runs. 



Figure A.8.1. — Starting longitude 0° results. 



Figure A.8.2. — Starting longitude 90° results. 
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A.9. Initial Covariance Case 9 Results 


Figures A.9.1 and A. 9.2 show the range covariance results for the 0° and 90°E starting longitude 
conditions, respectively, when the initial covariance is operating in Case 9 parameters (Cartesian position 
components =10 nr, Cartesian velocity components = 10 (m/s) 2 ). Blue lines represent range covariance 
results for the current system based on pseudo-range and ADR measurements from the 6 MS ground 
stations. Red lines represent range covariance results for MSysl, combining the current system with laser 
ranging measurements from the 6 MS stations and 2 additional ground stations. Similarly, black lines 
represent range covariance results for MSys2; green lines represent range covariance results for MSys3; 
cyan lines represent range covariance results for MSys4; yellow lines represent range covariance results 
for MSys5. Note that there are 10 individual lines for each color, derived from the 10 noise profile runs. 
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Figure A.9.1. — Starting longitude 0° results. 



Figure A. 9. 2. — Starting longitude 90° results. 
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